suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(tinytable))
suppressPackageStartupMessages(library(rairtable))
suppressPackageStartupMessages(library(lme4))
suppressPackageStartupMessages(library(broom))
sample_metadata <- read_tsv(paste0("data/metadata.tsv")) %>%
mutate(Extraction=factor(Extraction, levels=c("ZYMO","DREX","EHEX"))) %>%
rename(dataset=Dataset)
extraction_data <- airtable("tblBcTZcRG1E9wsGO", "appQpr6MxnaiVHsHy") %>% #get base ID from Airtable browser URL
read_airtable(., fields = c("EX DNA ng","Datasets_flat"), id_to_col = TRUE) %>% #get 3 columns from MAGs table
rename(id=1,extract=2,dataset=3) %>%
filter(dataset %in% sample_metadata$dataset) %>%
select(dataset,extract)
library_data <- airtable("tblo6AuYpxbbGw9gh", "appQpr6MxnaiVHsHy") %>% #get base ID from Airtable browser URL
read_airtable(., fields = c("Required PCR cycles","Datasets_flat"), id_to_col = TRUE) %>% #get 3 columns from MAGs table
rename(id=1,pcr=2,dataset=3) %>%
filter(dataset %in% sample_metadata$dataset) %>%
select(dataset,pcr)
indexing_data <- airtable("tblhfsiR4NI9XJQG0", "appQpr6MxnaiVHsHy") %>% #get base ID from Airtable browser URL
read_airtable(., fields = c("Adaptors (nM)","Library (nM)","Datasets_flat"), id_to_col = TRUE) %>% #get 3 columns from MAGs table
rename(id=1,adaptors=2,library=3,dataset=4) %>%
filter(dataset %in% sample_metadata$dataset) %>%
mutate(ratio=library/(adaptors+library)) %>%
select(dataset,adaptors,library,ratio)
preprocessing_data <- airtable("tblJfLRU2FIVz37Y1", "appQpr6MxnaiVHsHy") %>% #get base ID from Airtable browser URL
read_airtable(., fields = c("bases_pre_fastp","bases_post_fastp","host_bases","metagenomic_bases","singlem_fraction","C","EHI_plaintext"), id_to_col = TRUE) %>% #get 3 columns from MAGs table
rename(id=1,nonpareil=C,dataset=EHI_plaintext,microbial_fraction=singlem_fraction) %>%
filter(dataset %in% sample_metadata$dataset) %>%
select(dataset,bases_pre_fastp,bases_post_fastp,bases_post_fastp,metagenomic_bases,microbial_fraction,nonpareil)
assembly_data <- airtable("tblG6ZIvkYN844I97", "appQpr6MxnaiVHsHy") %>% #get base ID from Airtable browser URL
read_airtable(., fields = c("EHI_number_api","assembly_length","N50","L50","num_contigs","num_bins","metagenomic_bases"), id_to_col = TRUE) %>% #get 3 columns from MAGs table
rename(id=1,assembly_mapped_bases=metagenomic_bases,dataset=EHI_number_api) %>%
filter(dataset %in% sample_metadata$dataset) %>%
select(dataset,assembly_length,N50,L50,num_contigs,num_bins,assembly_mapped_bases)
mapping_data <- airtable("tblWDyQmM9rQ9wq57", "appWbHBNLE6iAsMRV") %>% #get base ID from Airtable browser URL
read_airtable(., fields = c("EHI_sample_static","MAG_mapping_percentage"), id_to_col = TRUE) %>% #get 3 columns from MAGs table
rename(id=1,dataset=EHI_sample_static) %>%
filter(dataset %in% sample_metadata$dataset) %>%
select(dataset,MAG_mapping_percentage)
all_data <- sample_metadata %>%
mutate(Taxon=factor(Taxon,levels=c("Amphibian","Bird","Mammal","Reptile","Control"))) %>%
left_join(extraction_data, by=join_by(dataset==dataset)) %>%
left_join(library_data, by=join_by(dataset==dataset)) %>%
left_join(indexing_data, by=join_by(dataset==dataset)) %>%
left_join(preprocessing_data, by=join_by(dataset==dataset)) %>%
left_join(assembly_data, by=join_by(dataset==dataset)) %>%
left_join(mapping_data, by=join_by(dataset==dataset))
Total amount of DNA extracted from the 150 ul subset of the bead-beaten sample.
all_data %>%
group_by(Taxon,Extraction) %>%
summarise(value = sprintf("%.0f±%.0f", mean(extract), sd(extract))) %>%
pivot_wider(names_from = Extraction, values_from = value) %>%
tt(caption = "Mean and standard deviation of total DNA nanograms")
| Taxon | ZYMO | DREX | EHEX |
|---|---|---|---|
| Amphibian | 297±350 | 340±475 | 1138±1195 |
| Bird | 71±134 | 11±11 | 44±45 |
| Mammal | 448±432 | 256±224 | 867±730 |
| Reptile | 102±71 | 162±134 | 250±179 |
| Control | 0±0 | 2±3 | 1±0 |
all_data %>%
ggplot(aes(x=Extraction,y=extract))+
geom_boxplot() +
facet_grid(. ~ Taxon, scales = "free") +
labs(y="DNA yield (ng)",x="Extraction method")
all_data %>%
filter(Taxon != "Control") %>%
lmerTest::lmer(extract ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
broom.mixed::tidy() %>%
tt()
| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 229.412625000 | 129.0917 | 1.7771296 | 18.62337 | 0.091883533 |
| fixed | NA | ExtractionDREX | -37.448812500 | 100.5929 | -0.3722808 | 59.99998 | 0.710995523 |
| fixed | NA | ExtractionEHEX | 345.333000000 | 100.5929 | 3.4329752 | 59.99998 | 0.001087593 |
| ran_pars | Sample | sd__(Intercept) | 0.003841839 | NA | NA | NA | NA |
| ran_pars | Species | sd__(Intercept) | 373.178642666 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 348.464099979 | NA | NA | NA | NA |
Number of PCR cycles required for reaching the plateau phase of the indexing PCR.
all_data %>%
group_by(Taxon,Extraction) %>%
summarise(value = sprintf("%.1f±%.1f", mean(pcr), sd(pcr))) %>%
pivot_wider(names_from = Extraction, values_from = value) %>%
tt(caption = "Mean and standard deviation of optimal number of PCR cycles")
| Taxon | ZYMO | DREX | EHEX |
|---|---|---|---|
| Amphibian | 14.2±2.1 | 10.8±3.3 | 10.2±2.6 |
| Bird | 18.0±2.5 | 18.5±3.9 | 14.8±2.6 |
| Mammal | 11.0±4.0 | 10.2±1.9 | 10.3±2.2 |
| Reptile | 11.7±4.2 | 10.3±3.7 | 12.0±3.5 |
| Control | 20.0±2.8 | 19.8±0.5 | 22.0±4.2 |
all_data %>%
ggplot(aes(x=Extraction,y=pcr))+
geom_boxplot() +
facet_grid(. ~ Taxon, scales = "free") +
labs(y="Optimal number of PCR cycles",x="Extraction method")
all_data %>%
filter(Taxon != "Control") %>%
lmerTest::lmer(pcr ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
broom.mixed::tidy() %>%
tt()
| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 13.708333 | 0.9696194 | 14.137850 | 22.50455 | 1.117002e-12 |
| fixed | NA | ExtractionDREX | -1.250000 | 0.8896915 | -1.404981 | 60.00000 | 1.651832e-01 |
| fixed | NA | ExtractionEHEX | -1.875000 | 0.8896915 | -2.107472 | 60.00000 | 3.926431e-02 |
| ran_pars | Sample | sd__(Intercept) | 0.000000 | NA | NA | NA | NA |
| ran_pars | Species | sd__(Intercept) | 2.555902 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 3.081982 | NA | NA | NA | NA |
#Total GB
all_data %>% summarise(sum(bases_pre_fastp)) / 1000000000
## sum(bases_pre_fastp)
## 1 429.0952
#Total reads
all_data %>% summarise(sum(bases_pre_fastp)) / 300
## sum(bases_pre_fastp)
## 1 1430317433
#Average GBs per sample
all_data %>% summarise(mean(bases_pre_fastp)) / 1000000000
## mean(bases_pre_fastp)
## 1 5.36369
all_data %>% summarise(sd(bases_pre_fastp)) / 1000000000
## sd(bases_pre_fastp)
## 1 2.703487
#Average reads per sample
all_data %>% summarise(mean(bases_pre_fastp)) / 300
## mean(bases_pre_fastp)
## 1 17878968
all_data %>% summarise(sd(bases_pre_fastp)) / 300
## sd(bases_pre_fastp)
## 1 9011624
all_data %>%
group_by(Taxon,Extraction) %>%
summarise(value = sprintf("%.1f±%.1f", mean(bases_post_fastp / 1000000000), sd(bases_post_fastp / 1000000000))) %>%
pivot_wider(names_from = Extraction, values_from = value) %>%
tt(caption = "Mean and standard deviation of sequencing depth (GB)")
| Taxon | ZYMO | DREX | EHEX |
|---|---|---|---|
| Amphibian | 4.0±1.5 | 4.2±2.3 | 4.7±0.4 |
| Bird | 3.7±1.9 | 3.4±2.5 | 2.7±1.9 |
| Mammal | 5.5±3.3 | 4.6±2.2 | 3.8±2.5 |
| Reptile | 6.1±2.3 | 5.7±1.4 | 5.0±1.9 |
| Control | 0.0±0.0 | 0.5±0.6 | 2.1±2.7 |
all_data %>%
ggplot(aes(x=Extraction,y=bases_pre_fastp))+
geom_boxplot() +
facet_grid(. ~ Taxon, scales = "free") +
labs(y="DNA yield (ng)",x="Extraction method")
all_data %>%
filter(Taxon != "Control") %>%
lmerTest::lmer(bases_post_fastp ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
broom.mixed::tidy() %>%
tt()
| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 4844471692 | 451697359 | 10.7250388 | 3.202405e+01 | 3.943247e-12 |
| fixed | NA | ExtractionDREX | -386734817 | 487241483 | -0.7937231 | 3.392332e+20 | 4.273567e-01 |
| fixed | NA | ExtractionEHEX | -786257695 | 487241483 | -1.6136920 | 8.048610e+20 | 1.065942e-01 |
| ran_pars | Sample | sd__(Intercept) | 1195869523 | NA | NA | NA | NA |
| ran_pars | Species | sd__(Intercept) | 555777398 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 1687854008 | NA | NA | NA | NA |
all_data %>%
mutate(qf_bases=bases_post_fastp/bases_pre_fastp*100) %>%
group_by(Taxon,Extraction) %>%
summarise(value = sprintf("%.1f±%.1f", mean(qf_bases), sd(qf_bases))) %>%
pivot_wider(names_from = Extraction, values_from = value) %>%
tt(caption = "Mean and standard deviation of quality-filtered proportion of reads")
| Taxon | ZYMO | DREX | EHEX |
|---|---|---|---|
| Amphibian | 84.6±1.5 | 91.1±4.8 | 87.5±4.0 |
| Bird | 70.6±16.5 | 62.5±29.3 | 66.7±19.3 |
| Mammal | 92.2±2.4 | 89.0±5.1 | 91.3±1.8 |
| Reptile | 89.9±6.6 | 90.5±7.4 | 88.3±7.6 |
| Control | 3.3±2.3 | 9.8±11.5 | 27.5±3.4 |
all_data %>%
mutate(qf_bases=bases_post_fastp/bases_pre_fastp*100) %>%
ggplot(aes(x=Extraction,y=qf_bases))+
geom_boxplot() +
facet_grid(. ~ Taxon, scales = "free") +
labs(y="DNA yield (ng)",x="Extraction method")
all_data %>%
mutate(qf_bases=bases_post_fastp/bases_pre_fastp*100) %>%
filter(Taxon != "Control") %>%
lmerTest::lmer(qf_bases ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
broom.mixed::tidy() %>%
tt()
| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 84.3251792 | 3.587497 | 23.5052981 | 18.58420 | 2.789261e-15 |
| fixed | NA | ExtractionDREX | -1.0626469 | 2.798710 | -0.3796917 | 48.00002 | 7.058490e-01 |
| fixed | NA | ExtractionEHEX | -0.8627821 | 2.798710 | -0.3082785 | 48.00002 | 7.592043e-01 |
| ran_pars | Sample | sd__(Intercept) | 6.6948031 | NA | NA | NA | NA |
| ran_pars | Species | sd__(Intercept) | 9.2214285 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 9.6950144 | NA | NA | NA | NA |
host_stats <- read_tsv("data/host_stats.tsv")
host_stats %>%
left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
group_by(Taxon,Extraction) %>%
summarise(value = sprintf("%.1f±%.1f", mean(duplicates), sd(duplicates))) %>%
pivot_wider(names_from = Extraction, values_from = value) %>%
tt(caption = "Mean and standard deviation of fraction of duplicated reads")
| Taxon | ZYMO | DREX | EHEX |
|---|---|---|---|
| Amphibian | 0.3±0.2 | 0.2±0.2 | 0.2±0.2 |
| Bird | 0.8±0.3 | 0.9±0.1 | 0.6±0.4 |
| Control | 1.0±0.0 | 0.9±0.0 | 1.0±0.0 |
| Mammal | 0.4±0.4 | 0.2±0.1 | 0.2±0.2 |
| Reptile | 0.5±0.4 | 0.3±0.3 | 0.4±0.4 |
host_stats %>%
left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
ggplot(aes(x=Extraction,y=duplicates))+
geom_boxplot() +
facet_grid(. ~ Taxon, scales = "free") +
labs(y="Duplication rate",x="Extraction method")
host_stats %>%
left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
filter(Taxon != "Control") %>%
lmerTest::lmer(duplicates ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
broom.mixed::tidy() %>%
tt()
| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 0.48796262 | 0.08012126 | 6.090302 | 20.71993 | 5.105684e-06 |
| fixed | NA | ExtractionDREX | -0.08741146 | 0.06897524 | -1.267287 | 60.00000 | 2.099492e-01 |
| fixed | NA | ExtractionEHEX | -0.15230321 | 0.06897524 | -2.208085 | 60.00000 | 3.107140e-02 |
| ran_pars | Sample | sd__(Intercept) | 0.00000000 | NA | NA | NA | NA |
| ran_pars | Species | sd__(Intercept) | 0.22019874 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 0.23893725 | NA | NA | NA | NA |
host_stats %>%
left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
group_by(Taxon,Extraction) %>%
summarise(value = sprintf("%.1f±%.1f", mean(coverage_depth), sd(coverage_depth))) %>%
pivot_wider(names_from = Extraction, values_from = value) %>%
tt(caption = "Mean and standard deviation of fraction of duplicated reads")
| Taxon | ZYMO | DREX | EHEX |
|---|---|---|---|
| Amphibian | 0.0±0.0 | 0.0±0.0 | 0.0±0.0 |
| Bird | 0.4±0.8 | 0.6±0.5 | 0.8±0.8 |
| Control | 0.0±0.0 | 0.0±0.0 | 0.0±0.0 |
| Mammal | 0.7±1.2 | 0.3±0.4 | 0.5±1.0 |
| Reptile | 0.2±0.3 | 0.1±0.1 | 0.1±0.1 |
host_stats %>%
left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
ggplot(aes(x=Extraction,y=coverage_depth))+
geom_boxplot() +
facet_grid(. ~ Taxon, scales = "free") +
labs(y="Depth of coverage",x="Extraction method")
host_stats %>%
left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
filter(Taxon != "Control") %>%
lmerTest::lmer(coverage_depth ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
broom.mixed::tidy() %>%
tt()
| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 0.329125000 | 0.1331063 | 2.47264850 | 20.58494 | 0.02222872 |
| fixed | NA | ExtractionDREX | -0.089791667 | 0.1144640 | -0.78445305 | 48.00000 | 0.43662873 |
| fixed | NA | ExtractionEHEX | -0.002791667 | 0.1144640 | -0.02438903 | 48.00000 | 0.98064341 |
| ran_pars | Sample | sd__(Intercept) | 0.408634577 | NA | NA | NA | NA |
| ran_pars | Species | sd__(Intercept) | 0.224731208 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 0.396515070 | NA | NA | NA | NA |
host_stats %>%
left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
group_by(Taxon,Extraction) %>%
summarise(value = sprintf("%.1f±%.1f", mean(coverage_breadth), sd(coverage_breadth))) %>%
pivot_wider(names_from = Extraction, values_from = value) %>%
tt(caption = "Mean and standard deviation of depth of host genome coverage")
| Taxon | ZYMO | DREX | EHEX |
|---|---|---|---|
| Amphibian | 0.0±0.0 | 0.0±0.0 | 0.0±0.0 |
| Bird | 0.6±0.5 | 3.2±4.4 | 8.9±13.9 |
| Control | 0.0±0.0 | 0.0±0.0 | 0.0±0.0 |
| Mammal | 5.7±5.9 | 10.2±16.4 | 15.1±26.4 |
| Reptile | 4.9±7.5 | 3.0±5.4 | 2.9±3.7 |
host_stats %>%
left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
ggplot(aes(x=Extraction,y=coverage_breadth))+
geom_boxplot() +
facet_grid(. ~ Taxon, scales = "free") +
labs(y="Breadth of coverage",x="Extraction method")
host_stats %>%
left_join(sample_metadata,by=join_by(dataset==dataset)) %>%
filter(Taxon != "Control") %>%
lmerTest::lmer(coverage_breadth ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
broom.mixed::tidy() %>%
tt()
| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 2.799167 | 2.252826 | 1.2425133 | 20.82558 | 0.22785670 |
| fixed | NA | ExtractionDREX | 1.301250 | 1.956426 | 0.6651158 | 48.00000 | 0.50915975 |
| fixed | NA | ExtractionEHEX | 3.918750 | 1.956426 | 2.0030144 | 48.00000 | 0.05084021 |
| ran_pars | Sample | sd__(Intercept) | 7.118462 | NA | NA | NA | NA |
| ran_pars | Species | sd__(Intercept) | 3.549767 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 6.777259 | NA | NA | NA | NA |
Nonpareil estimate of the metagenomic complexity after removing host DNA.
all_data %>%
group_by(Taxon,Extraction) %>%
summarise(value = sprintf("%.1f±%.1f", mean(nonpareil), sd(nonpareil))) %>%
pivot_wider(names_from = Extraction, values_from = value) %>%
tt(caption = "Mean and standard deviation of breadth of host genome coverage")
| Taxon | ZYMO | DREX | EHEX |
|---|---|---|---|
| Amphibian | 0.9±0.1 | 0.8±0.1 | 0.8±0.1 |
| Bird | 0.9±0.1 | 1.0±0.0 | 0.8±0.4 |
| Mammal | 0.8±0.2 | 0.8±0.1 | 0.7±0.3 |
| Reptile | 0.9±0.1 | 0.9±0.0 | 0.9±0.1 |
| Control | 0.0±0.0 | 0.5±0.6 | 0.5±0.7 |
all_data %>%
ggplot(aes(x=Extraction,y=nonpareil))+
geom_boxplot() +
facet_grid(. ~ Taxon, scales = "free") +
labs(y="Nonpareil completeness",x="Extraction method")
all_data %>%
filter(Taxon != "Control") %>%
lmerTest::lmer(nonpareil ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
broom.mixed::tidy() %>%
tt()
| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 0.885552853 | 0.03450229 | 25.6664948 | 40.89542 | 7.753071e-27 |
| fixed | NA | ExtractionDREX | 0.005536892 | 0.04336611 | 0.1276779 | 48.00020 | 8.989373e-01 |
| fixed | NA | ExtractionEHEX | -0.112193866 | 0.04336611 | -2.5871326 | 48.00020 | 1.276294e-02 |
| ran_pars | Sample | sd__(Intercept) | 0.059565487 | NA | NA | NA | NA |
| ran_pars | Species | sd__(Intercept) | 0.035030813 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 0.150224598 | NA | NA | NA | NA |
Variance partitioning metrics are derived from community_analysis.Rmd.
alpha_data <- list.files(path = "results", pattern = "^alpha_.*\\.tsv$", full.names = TRUE) %>%
map_df(~ read_tsv(.)) %>%
left_join(sample_metadata,by= join_by(dataset==dataset))
alpha_data %>%
pivot_longer(!c(dataset,Library,Species,Taxon,Sample,Extraction), names_to = "metric", values_to = "value") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
ggplot(aes(x=Extraction,y=value))+
geom_boxplot() +
facet_grid(metric ~ Taxon, scales = "free")
alpha_data %>%
lmerTest::lmer(richness ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
broom.mixed::tidy() %>%
tt()
| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 58.43750 | 12.008400 | 4.8663854 | 8.87726 | 0.0009237868 |
| fixed | NA | ExtractionDREX | 2.18750 | 4.698926 | 0.4655319 | 32.00001 | 0.6447035294 |
| fixed | NA | ExtractionEHEX | 2.43750 | 4.698926 | 0.5187355 | 32.00001 | 0.6075141409 |
| ran_pars | Sample | sd__(Intercept) | 15.62413 | NA | NA | NA | NA |
| ran_pars | Species | sd__(Intercept) | 30.71216 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 13.29057 | NA | NA | NA | NA |
alpha_data %>%
lmerTest::lmer(neutral ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
broom.mixed::tidy() %>%
tt()
| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 29.71528470 | 5.751141 | 5.16685008 | 8.745982 | 0.0006458923 |
| fixed | NA | ExtractionDREX | -1.59472963 | 2.085900 | -0.76452835 | 31.999998 | 0.4501538272 |
| fixed | NA | ExtractionEHEX | -0.04935512 | 2.085900 | -0.02366131 | 31.999998 | 0.9812697035 |
| ran_pars | Sample | sd__(Intercept) | 9.00571703 | NA | NA | NA | NA |
| ran_pars | Species | sd__(Intercept) | 14.37531297 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 5.89981588 | NA | NA | NA | NA |
alpha_data %>%
lmerTest::lmer(phylogenetic ~ Extraction + (1 | Sample) + (1 | Species), data = ., REML = FALSE) %>%
broom.mixed::tidy() %>%
tt()
| effect | group | term | estimate | std.error | statistic | df | p.value |
|---|---|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 4.5870074 | 0.4577012 | 10.021838 | 8.589824 | 4.977993e-06 |
| fixed | NA | ExtractionDREX | 0.2122362 | 0.1485377 | 1.428837 | 32.000064 | 1.627406e-01 |
| fixed | NA | ExtractionEHEX | 0.2230529 | 0.1485377 | 1.501658 | 32.000064 | 1.429897e-01 |
| ran_pars | Sample | sd__(Intercept) | 1.1278928 | NA | NA | NA | NA |
| ran_pars | Species | sd__(Intercept) | 0.9754991 | NA | NA | NA | NA |
| ran_pars | Residual | sd__Observation | 0.4201282 | NA | NA | NA | NA |
Variance partitioning metrics are derived from community_analysis.Rmd.
variance_data <- list.files(path = "results", pattern = "^var_.*\\.tsv$", full.names = TRUE) %>%
map_df(~ read_tsv(.))
variance_data %>% summarise(mean=mean(r2),sd=sd(r2))
## # A tibble: 1 × 2
## mean sd
## <dbl> <dbl>
## 1 0.102 0.0768
variance_data %>%
left_join(sample_metadata %>% select(Species,Taxon) %>% unique(),by=join_by(species==Species)) %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic"))) %>%
ggplot(aes(x=Taxon,y=r2)) +
geom_boxplot()+
ylim(0,1)+
facet_grid(. ~ metric, scales = "free")
variance_data %>%
group_by(metric) %>%
summarise(mean=mean(r2),sd=sd(r2)) %>%
tt()
| metric | mean | sd |
|---|---|---|
| neutral | 0.06201935 | 0.04744380 |
| phylogenetic | 0.07381277 | 0.06884687 |
| richness | 0.16970818 | 0.06626388 |